High Pressure Thermoelasticity of Body-centered Cubic Tantalum 
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We have investigated the thermoelasticity of body-centered cubic (bcc) tantalum from first prin- 
ciples by using the linearized augmented plane wave (LAPW) and mixed-basis pseudopotential 
methods for pressures up to 400 GPa and temperatures up to 10000 K. Electronic excitation contri- 
butions to the free energy were included from the band structures, and phonon contributions were 
included using the particle-in-a-cell (PIC) model. The computed elastic constants agree well with 
available ultrasonic and diamond anvil cell data at low pressures, and shock data at high pressures. 
The shear modulus C44 and the anisotropy change behavior with increasing pressure around 150 GPa 
because of an electronic topological transition. We find that the main contribution of temperature 
to the elastic constants is from the thermal expansivity. The PIC model in conjunction with fast 
self-consistent techniques is shown to be a tractable approach to studying thermoelasticity. 

PACS numbers: 62.20.Dc, 62.20.-x, 71.20.Be, 46.25.Hf, 46.25.-y 



Single crystal elastic constants of solids at high pres- 
sures and temperatures are essential in order to predict 
and understand material response, strength, mechani- 
cal stability, and phase transitions. We have studied 
the high pressure and temperature elastic constants of 
body-centered cubic (bcc) tantalum, a group V tran- 
sition metal, from first principles. Because of its high 
structural mechanical, thermal and .chemical stability, Ta 
is a useful high pressure standard.El Ta has a very high 
melting temperature, 3269 K at ambient pressure. Bcc 
Ta is stable to 174 GPa, according to diamondUanvil-cell 
experiments.^ Shock compression experimentsa show no 
transition other than melting (at around 300 GPa). Its 
stability makes Ta an ideal material for understanding 
the generic behavior of transition metals under compres- 
sion, without the complication of phase transitions. Re- 
cently, its static properties were studied by full-potential 
LMTO calculationsB and the thermal equation of state 
was reported.u 

The three elastic constants, en, cyi and C44, completely 
describe the elastic behaviour of a cubic crystal. A more 
convenient set for computations are C44 and two linear 
combinations, K and c s . The bulk modulus, 

K = ( Cll + 2c 12 )/3, (1) 

is the resistance to deformation by an uniform hydro- 
static pressure; the shear constant, 

c s = (c n - c 12 )/2, (2) 

is the resistance to shear deformation across the (110) 
plane in the [110] direction, and C44 is the resistance to 
shear deformation across the (100) plane in the [010] di- 
rection. The bulk-modulus K was determinedjirom the 
equation of statejj using the Vinet equation.tffi We ob- 
tained the shear moduli by straining the bcc lattice at 
fixed volumes using volume conserving tetragonal and or- 



thorhombic strains for c s and C44 respectively, and com- 
puting the free energy as a function of strain. c s was 
obtained by applying the following isochoric strain 
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where <5 is the magnitude of the strain. Then the strain 
energy is 

F(S) = F(0) + 6c s V5 2 + O(S 3 ), (4) 

where F(0) is the free energy of the unstrained system 
and V is its volume. Similarly, C44 was calculated from 
the following strain 

{ 6 \ 
e = { 5 (5) 
\ 5 2 /(l-5 2 ) J 

with the corresponding strain energy 

F(5) = F(0) + 2 Cii V6 2 + 0{S i ). (6) 

The quadratic coefficients of strain energy gives the elas- 
tic constants. First order terms due to the initial stress 
(hydrostatic pressure)Q were eliminated by applying iso- 
choric strains. Then, the elastic constants cn and C12 
were obtained from c s and K . 

We assume that the HeLmholtz free energy of the sys- 
tem can be separated asEnLTJ: 

F(V, T) = E Q (V) + F el (V, T) + F vib (V, T) (7) 

where E$(V) is the static zero temperature energy, 
F e i(V,T) is the electronic contribution, and F v u,(V,T) 
is the vibrational contribution to the free energy. Our 
computational procedure is based on density functional 
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theory (DFT) generalized to finite temperatures by the 
Mermin theorem.Ea The charge density is temperature 
dependent through occupation numbers according to the 
Fermi-Dirac distribution, giving the electronic entropy 
from 

Sei = Y, & ln fi + i 1 ~ /<)Mi - h) (8) 

where fa = fi(E — Ep,T)) is the Fermi occupation at T 
for each state i. The variations of fi with temperature 
were included from the self-consistent band structures 
calculated at an electronic temperature of 2000 K varying 
according to the Fermi-Dirac distribution. 

The electronic excitations, both the static energy and 
the electronic contribution to the free energy, were com- 
puted by using the full potential, linearized augmented 
plane wave (LAPW) method.EHtia The 5p,4/,5d and 6s 
states were treated as band states, and the deeper states 
were treated as soft core electrons. The generalized gra- 
dient approximation (GGA)ta was used for the exchange- 
correlation potential. The convergence of strain energies 
with respect to the Brillouin zone integration was care- 
fully checked by repeating the calculations for 16x16x16 
and 24x24x24 meshes at V=16.82 A 3 and we found at 
most 2 GPa (3 %) difference both for c s aneLc44. Hence, 
we used 16x16x16 special k-points meshestj in the full 
Brillouin zone giving 344 and 612 k-points within IBZ 
of tetragonal and orthorhombic lattice respectively. The 
convergence parameter RK rnax was 9 giving about 1800 
planewaves and 200 basis functions per atom at zero pres- 
sure. 

The vibrational free-energy was obtained within the 
particle-in-a-cell model (PIC)Ej by using an accu=. 
rate pseudopotential mixed-basis total energy methodhll 
which is computationally more efficient than the LAPW 
calculations. In PIC, an atom is displaced in its Wigner- 
Seitz cell in the potential field of all the other atoms fixed 
at their equilibrium positions, i.e. the ideal, static lattice 
except for the wanderer atom. The partition function, 
and hence the free energy is calculated from this poten- 
tial energy surface via an integral over the position of a 
single atom inside the Wigner-Seitz cell. The PIC model 
is essentially an anharmonic Einstein model, and the 3N 
dimensional partition function is reduced to a simple 3D 
integral.ua The advantage of the cell model over lattice 
dynamics based on the quasiharmonic approximation is 
that anharmonic contributions from the potential-energy 
of the system have been included exactly without a per- 
turbation expansion. On the other hand, since we used 
the classical partition function, and the interatomic cor- 
relations between the motions of different atoms is ig- 
nored, it is only valid at high temperatures above the 
Debye temperature (245 K in Ta.Efh-Since the vacancy 
formation energy is very high in Ta,tj spontaneous for- 
mation of defects is only important after the melting tem- 
perature. 
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FIG. 1. Static elastic constants of Ta as a function of 
pressure. Open-Bewares are ultrasonic experimental data of 
Katahara et a£J'E3, and the dotted lines show the-initial 
slopes. Open symbols are SAX data of Cynn and Yoocj. The 
anisotropy ratio is shown in the inset. 

For the PIC computations, a supercell with 54 atoms 
was used. The pseudopotential mixed-basis calcula- 
tions were, carried out on this 54 atoms supercell us- 
ing LDAEJ for exchange-correlations effects and 2x2x2 
k-point mesh resulting 4 special k points for BZ in- 
tegrations. A semi-relatiyiS|tic, nonlocal and norm- 
conserving Troullier-MartinsO pseudopotential (with as- 
sociated pseudo-atomic orbitals) with non-linear core 
correctional was used for the Ta atoms as described in 
detail ip our previous study of thermal equation of state 
of Ta.u After checking the energy convergence, 550 eV 
and 60 eV are used as planewave energy cutoffs for the 
expansion of the pseudo-atomic orbitals as well as FFT 
grid and low-energy plane waves for additional degrees 
of freedom in basis set respectively. The canonical par- 
tition function was computed from the potential energy 
surface as a function of displacements, of wanderer atom 
along special symmetry directions. 3o We used 2 and 4 
special directions for tetragonal and orthorhombic distor- 
tions, respectivefv—which integrates exactly up to I = 6 
lattice harmonicsrJ The potential energy was calculated 
at 4-6 different displacements along each of these special 
directions, and was fit to an even polynomials up to or- 
der six. Details of the. all computational parameters were 
described previously.u 

The static elastic constants as functions of pressure are 
presented in Fig. [I] and Table El The zero pressure values 
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and initial slopes are in good agreement wittuthe ultra- 
sonic experimental data of Katahara et aZEj'Ej Similarly, 
comparison with recent SAX (stoss/angle-resolved x-ray 
diffraction) experimental datao up to 105 GPa shows 
good agreement for en and C44. Likewise, c\i agrees well 
at low pressures, but deviates with increasing pressure. 
This may be due to the assumed isostress condition for 
experimental data analysis for all pressures, or due to the 
large uncertainty on measured deviatoric stress at high 
pressures. Note that, the initial slope of ultrasonic data 
agrees very well with our calculated C\2- The anisotropy 
ratio, A = c^jc s (inset Fig. |l|) first decreases from 1.57 
to 0.9 with increasing pressure, and then its slope re- 
verses and it increases with increasing pressure. This is 
due to the changes in C44. The change in behavior of C44 
around 150 GPa is due to the. electronic transition evi- 
dent in the equation of state.Q This indicates that elas- 
tic constants can be much more sensitive to changes in 
the Fermi surface than the equation of state, where the 
electronic transition was not apparent without examining 
small residuals in the equation of state fit. 

The elastic constants of Ta as functions of temperature 
at various pressures are presented in Fig. || and Table |J. 
In order to compare with experimental data, the com- 
puted isothermal elastic constants (cjl are converted to 
adiabatic constants (cf ) according t< 
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where = J^k a k c Jki a k is the linear thermal expansion 
tensor, Cy is the specific heat and p is the density. For 
cubic crystals, eq. O simplifies to 



FIG. 2. The elastic constants of bcc Ta as a function of 
temperature at different pressure from GPa (lowest curve) 
to 400 GPa (uppermost curve) with 50 GPa interval, a) shear 
modulus c s , b) shear modulus C44 and c) adiabatic bulk mod- 
ulus i££p. Dotted lines are the experimental data from Walker 
etalM 
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1965.52 


1181.94 



TABLE I. The static elastic constants for bcc Tantalum. 
All elastic constants as well as pressure units are GPa. 



with a is the thermal expansion coefficient, 7 is the 
Griineisen parameter, and Kt is the isothermal bulk 
modulus. The thermodynamic parameters were com- 
puted, self-consistently from the thermal equation of 
state.Q The correction is zero for C44 and c s . A increases 
with temperature but decreases with pressure; at 3000 K 
it decreases from 5 % to 1 % for pressures 50 GPa to 
400 GPa for bulk modulus K T , and at 10000 K A is 
29 % and 3 % for same pressures. 



T (K) 



C44 



,.T 
Cll 



Cll 



Cl2 



cf 2 



K T K s 



44.05 70.26 249.68 249.68 161.59 161.59 190.95 190.95 

947 43.58 56.24 221.71 233.41 134.55 146.27 163.61 175.33 

2053 38.61 59.39 189.66 211.95 112.44 134.72 138.18 160.46 

3000 31.57 62.48 162.55 192.74 99.41 129.60 120.46 150.65 

3947 28.72 64.47 142.15 179.77 84.70 122.32 103.85 141.47 

5052 18.34 59.70 107.55 154.15 70.87 117.47 83.10 129.70 

6000 0.31 62.05 64.91 120.03 64.29 119.40 64.50 119.61 

TABLE II. The elastic constants for bcc Tantalum at var- 
ious temperatures. All elastic constants units are GPa. 
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The shear moduli c s and C44 and adiabatic bulk modtu. 
lus Ks agree well with the ultrasonic experimental dataE3 
up to 3000 K at zero pressure (Fig. |J). We find that all 
three moduli are primarily functions of volume, and ther- 
mal effects at constant volume are quite small except at 
the highest pressures. There is some softening of c s with 
increasing temperature for all pressures. C44, shows a 
slight softening at the zero pressure with increasing tem- 
perature, but they are rather flat for other pressures ex- 
cept than very high pressures. The adiabatic bulk mod- 
ulus Ks also softens slightly with temperature at low 
pressures but becomes flat with increasing pressure. 
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FIG. 3. The anisotropy ratio of elastic constants of Ta as 
a function of temperature at different pressures from GPa 
to 400 GPa with 50 GPa interval. 

The anisotropy ratio, A = 04,4,/ c s , is presented as a 
function of temperature for various pressures at Fig. |^. 
A increases with increasing temperature at all pressures, 
but less drastically at high pressures. At lower pressures, 
this increase is divergent after certain temperature, since 
the softening of c s is large enough and it approaches 
zero. The reversal of the slope of A with pressure shifts 
to higher pressures with increasing temperature due to 
thermal expansivity, and occurs at a fixed volume. 

Sound velocities are related to the elastic constants by 
the Christoffel equatiorO 

{cijkinjUk - pv 2 S lj )u i = 0, (13) 

where Cijki is the elastic constants tensor, n is the prop- 
agation direction, u is the polarization vector and v is 
the velocity. Our elastic constants are those appropriate 
for the, equations of motion under hydrostatic reference 
stress. For [110] wave propagation direction in a cubic 
lattice, the longitudinal mode is 

pv 2 = (en +ci2 + 2c 44 )/2 (14) 

and two transverse modes are 

pv 2 = C44 (15) 

and 



pv 2 = (en - ci 2 )/2 = c s (16) 

polarized along [001] and [110] direction respectively. For 
polycrystalline sample, the average isotropic shear mod- 
ulus G can be determined from single crystal elastic con- 
stants according to the Voigt-Reuss-Hill schemeE], and 
the isotropically averaged aggregate velocities are given 
by 

v P = ((K + 4/3G)/pf 2 (17) 
v s = (Glpfl 2 (18) 
v B = {K/p) 1 ' 2 (19) 

where Vp, Vs, and vb are the compressional, shear and 
bulk sound velocities. The sound velocities of Ta along 
the Hugoniot calculated from elastic constants are shown 
in Fig. and are comparedr-with the shock sound veloc- 
ity data from Brown et alE3 As seen in Fig. ||, there 
is excellent agreement with shock data. The calculated 
compressional velocity vp agrees very well with exper- 
imental data up to 200 GPa, and then after 300 GPa 
the bulk velocity vb matches the data well. This is be- 
cause the shocked solid melts around 300 GPa, so the 
liquid velocity might be represented by Vb- The devia- 
tion between 200 GPa and 300 GPa is probably due to 
premelting effects. 
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FIG. 4. Sound velocities of Ta along the Hugoniot calcu- 
lated from elastic constants. Solid lines are the longitudinal 
and two transverse sound velocities in [110] direction from 
single crystal elastic constants. The polarization of the shear 
waves are given in brackets. The isotropic aggregate velocities 
are shown by dashed lines, vp, vb, and vs are the compres- 
sional, bulk, and shear sound velooities. Filled dots are the 
shock data from Brown and Shaneio. 

In conclusion, the elasticity of bec Ta is investigated 
from first principles for pressures up to 400 GPa and 
temperatures up to 10000 K. The calculated static elastic 
constants are in good agreement with available ultrasonic 
and SAX experimental data. The shear modulus C44 and 
the anisotropy ratio A change behaviour with increasing 
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pressure around 150 GPa. Although, the shear modulus 
c s softens with increasing temperature at all pressures, 
c 44 and Ks soften with temperature at low pressures but 
then they are rather flat at higher pressures. The main 
effect of temperature for the thermoelasticity of Ta is due 
to thermal expansivity. The calculated sound velocities 
along the Hugoniot shows an excellent agreement with 
shock-wave experimental data. 
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